function  [w,delta_w,u]=initial     %Set initial guess
global c a_t b p e t n_t n_p x_t y_t ar nth0 pwr;

c=1;                %-1/c*\delta u(x)+a_t(x) u(x)-f(x,u(x))=0.
pwr=3;              %nonlinear term in PDE f(x,u(x))=u(x)^pwr

a_t=zeros(n_t,1);          %-1/c*\delta u(x)+a_t(x) u(x)-f(x,u(x))=0.
%a_t=zeros(n_t,1);         %-1/c*\delta u(x)+a_t(x) u(x)-f(x,u(x))=0.

x_ctr=0;                 %-1. 0.25 2 for dumpbell
y_ctr=-0.85;
%delta_w=-ones(n_t,1);  %1 for convex down, -1 for convex up and 0 for no concern
delta_w=zeros(n_t,1);
for i=1:n_t
   dist=(((x_t(i)-x_ctr))^2+((y_t(i)-y_ctr))^2);
%  if dist <= 1.00001      
%  if x_t(i) <= 0             %east-west  |
%  if y_t(i) <= 0             %north-south -
%  if y_t(i) <= x_t(i)        %northwest-souteast /
%  if x_t(i)+y_t(i) <= 0      %northeast-southwest \
%  if (x_t(i) > 0 & y_t(i) > 0) | (x_t(i) < 0 & y_t(i) < 0)   % +
%   if (y_t(i) <= x_t(i) && x_t(i)+y_t(i) >= 0) | (x_t(i) < y_t(i) && x_t(i)+y_t(i) < 0)                         % X
   if(dist<=0.0225)
       delta_w(i)=cos(2*pi*dist/0.09);
   end
end

u=assempde(b,p,e,t,1,a_t',-delta_w'); %Solve w_p from Delta_w_t at vertices
w=pdeintrp(p,t,u);                    %w_t solution at triangular points
w=w';
delta_w=a_t.*w+delta_w;
